# This file contains the S-plus code that was used to implement the estimators and functionals in "A Measure of Asymmetry Based on a New Necessary and Sufficient Condition for Symmetry", Sankhya (2014)
# As a fully working example the code herein can reproduce Table 2 of the paper. 
# The same functions (with different input data) can be used in reproducing other figures of the same paper 
# More general use of the estimates is also feasible by supplying the corresponding data set.  Use of the functions is self - explanatory and is further facilitated by the comments at he preamble of each function.

remove(objects())
#
#module(envstats)
options(echo=FALSE)
sd2theta <- function(sd)
{
  return(sqrt(pi/2)/sd)
}

theta2sd <- function(theta)
{
  return(sqrt(pi/2)/theta)
}


dhalfnorm <- function(x, theta=sqrt(pi/2), log=FALSE)
{
  sd.norm = theta2sd(theta)
  
  if (log)
    d = ifelse(x<0, -Inf, log(2)+dnorm(x, mean=0, sd=sd.norm, log=TRUE))
  else
    d = ifelse(x<0, 0,  2*dnorm(x, mean=0, sd=sd.norm)) 

  return(d)
}

phalfnorm <- function(q, theta=sqrt(pi/2), lower.tail=TRUE, log.p=FALSE)
{
  sd.norm = theta2sd(theta)
  
  p = ifelse(q < 0, 0, 2*pnorm(q, mean=0, sd=sd.norm)-1)
  if (lower.tail == FALSE) p = 1-p
  if(log.p) p = log(p)
  
  return( p )
}

qhalfnorm <- function(p, theta=sqrt(pi/2), lower.tail=TRUE, log.p=FALSE)
{
  sd.norm = theta2sd(theta)
  
  if (log.p) p = exp(p)
  if (lower.tail == FALSE) p = 1-p
  q = ifelse(p < 0, NaN, qnorm((p+1)/2, mean=0, sd=sd.norm))
  
  return(q)
}

rhalfnorm <- function(n, theta=sqrt(pi/2))
{
  sd.norm = theta2sd(theta)
  return( abs(rnorm(n, mean=0, sd=sd.norm)) )
}
RSample<-function(s,dist, p1=0,p2=1)
{
	#cat("\n", dist, " ", p1, " ", p2, " ", s,"\n")
	switch(dist, weib = rweibull(s, shape=p1), 
					lognorm = rlnorm(s),
					norm = rnorm(s),
					uni = runif(s),
					cauchy = rcauchy(s, p1, p2)	,
					fnorm = abs(rnorm(s)	),
					normmixt = c(rnorm(ceiling(p1*s), 0, 1),rnorm(floor((1-p1)*s), 2,2))	#rnorm.mix(s, 0,1,2,2,p1)#			
					)

}

QSample<-function(s,dist, p1=0,p2=1)
{
	#cat("\n", dist)#, " ", p1, " ", p2, " ", s,"\n")
	switch(dist, weib = qweibull(s, shape=p1), 
					lognorm = qlnorm(s),
					norm = qnorm(s),
					uni = qunif(s),
					cauchy = qcauchy(s, p1, p2)	,
					fnorm = qhalfnorm(s),
					normmixt = quantile(c(rnorm(ceiling(p1*1000), 0, 1),rnorm(floor((1-p1)*1000), 2,2)), s)	#qnorm.mix(s, 0,1,2,2,p1)#			
					)

}

DSample<-function(s,dist, p1,p2)
{
	#cat("\n", dist, " ", p1, " ", p2, " ", s,"\n")
	switch(dist, weib = dweibull(s, shape=p1), 
					lognorm = dlnorm(s),
					norm = dnorm(s),
					uni = dunif(s),
					cauchy = dcauchy(s, p1, p2)	,
					fnorm = dhalfnorm(s),
					normmixt = p1 * dnorm(s,0,1)+(1-p1)*dnorm(s,2,2) #dnorm.mix(s, 0,1,2,2,p1)
					#ksmooth(s, kernel="normal", bandwidth=bandwidth.nrd(s), x.points=s[order(s)])$y + ksmooth(-s,  kernel="normal", bandwidth=bandwidth.nrd(s), x.points=s[order(s)])$y			
					)

}

DistSample<-function(s,dist, p1,p2)
{
	#cat("\n", dist, " ", p1, " ", p2, " ", s,"\n")
	switch(dist, weib = pweibull(s, shape=p1), 
					lognorm = plnorm(s),
					norm = pnorm(s),
					uni = punif(s),
					cauchy = pcauchy(s, p1, p2)	,
					fnorm = phalfnorm(s),
					normmixt = p1 * pnorm(s,0,1)+(1-p1)*pnorm(s,2,2) 
					#pnorm.mix(s, 0,1,2,2,p1)
					#ksmooth(s, kernel="normal", bandwidth=bandwidth.nrd(s), x.points=s[order(s)])$y + ksmooth(-s,  kernel="normal", bandwidth=bandwidth.nrd(s), x.points=s[order(s)])$y			
					)
}





Rho.p<-function(sample, p.param, dist, p1=0, p2=1)
{
	xi.p <- QSample(p.param, dist, p1,p2)
	lowerlim <- -Inf
	upperlim <- xi.p
	sampleuse<-sample[which(sample<=upperlim )]
	
	xpoints<-seq(min(sampleuse), max(sampleuse), length=50) #sampleuse[order(sampleuse)]
	pdfest<-    DSample(xpoints, dist, p1,p2)		
	cdfest<-DistSample(xpoints, dist, p1,p2)
	n<-length(xpoints)

	psi1<- sum(pdfest * cdfest)/n  # mean(pdfest * cdfest)
	psi2<-  sum(pdfest )/n  #mean(pdfest)
	psi3<-  sum(pdfest^2 )/n #mean(pdfest^2)
	
	denom<- p.param * psi3 - psi2^2
	rhop<- ifelse(denom<0, 0, 2*sqrt(3)/p.param * ( psi1 - p.param/2 * psi2 )/sqrt(p.param * psi3 - psi2^2))
		
	rhop
}

Rhostar.p<-function(sample, p.param, dist, p1, p2)
{
	xi.1minusp <- QSample(1-p.param, dist, p1,p2)
	lowerlim<- xi.1minusp	
	upperlim<- +Inf
	sampleuse<-sample[which( sample >= lowerlim)]
	xpoints<-seq(min(sampleuse), max(sampleuse), length=50)  #sampleuse[order(sampleuse)] #seq(min(sampleuse), max(sampleuse), length=50) #
	pdfest<-    DSample(xpoints, dist, p1,p2)		
	cdfest<-DistSample(xpoints, dist, p1,p2)
		
	n<-length(xpoints)

	psi1star<- sum(pdfest * cdfest)/n  # mean(pdfest * cdfest)
	psi2star<-  sum(pdfest )/n  #mean(pdfest)
	psi3star<-  sum(pdfest^2 )/n #mean(pdfest^2)
	
	denom<- p.param * psi3star - psi2star^2
	rhostarp<-ifelse(denom<0, 0, 2*sqrt(3)/p.param *  ( -psi2star + psi1star + p.param/2 * psi2star )/sqrt(p.param * psi3star - psi2star^2))

		
	
	rhostarp
}

SampleAsymCoef<-function(sample, p.param, dist, p1,p2)
{
	ro1<-Rho.p(sample, p.param, dist, p1, p2)
	ro2<-Rhostar.p(sample, p.param, dist, p1, p2)
	#abs(ro1+ro2)	
	
	ifelse1((ro1==0) || (ro2==0), 0, abs(ro1+ro2)	)
}



MaxAsymCoef<-function(SampleSize, dist, GridLength, p1=1, p2=0)
{
sample<- RSample( SampleSize, dist, p1, p2)
p<-seq(0.5, 0.99, length=GridLength)
r1<-sapply(1:GridLength, function(i, sample, p, dist, p1, p2) {SampleAsymCoef(sample, p[i], dist, p1,p2)}, sample, p, dist, p1,p2)
- 0.5 * sign(Rho.p(sample, 1, dist, p1, p2)) * max(r1)
}

Simulation<-function(ItTimes, SampleSize, dist, GridLength, p1=1, p2=0, TrueCoef)
#ItTimes: Number of iterations that the simulation will run
#SampleSize: Sample size to use
#dist: the distribution for which the simulation is performed
#GridLength: Over how many points between 0.5 and 0.99 the max will be calculated
#p1,p2: distribution parameters for all except normal mixture. For normal mixture use only p1: the proportion of mixture
{
	#Coefmatrix<-matrix(nrow=ItTimes, ncol=1)
	Coef<-sapply(1:ItTimes, function(i, SampleSize, dist, GridLength, p1, p2) {MaxAsymCoef(SampleSize, dist, GridLength, p1, p2)}, SampleSize, dist, GridLength, p1, p2)
#Coefmatrix<-Coef
result<-na.exclude(Coef)
#print.matrix(result, digits = 3)
c(mean((result-TrueCoef)^2), colMeans(result), stdev(result)^2)
}

SampleSize<-50
Iterations<-1000
set.seed(156)
#Matrix500All<-matrix(nrow=15, ncol=4)
#Matrix300All<-matrix(nrow=15, ncol=4)

Matrix300All<-matrix(nrow=15, ncol=4)


#colnames(Matrix500All)<-c("TrueCoef", "MSE","AverageValue", "variance") 
colnames(Matrix300All)<-c("TrueCoef", "MSE","AverageValue", "variance") 
#row.names(Matrix500All)<-c("W(2)","W(3)","W(4)","W(5)")#, "FN", "LN", "N", "NM1","NM2", "NM3","NM4", "NM5","NM6", "NM7","NM8")
row.names(Matrix300All)<-c("W(2)","W(3)","W(4)","W(5)", "FN", "LN", "N", "NM1","NM2", "NM3","NM4", "NM5","NM6", "NM7","NM8")
#Matrix500All[,1]<-c(0.42, 0.08, -0.07, -0.16)#, 0.984, 0.939, 0, 0.1,0.1,0.2,0.2,0.3,0.3,0.4,0.4)
Matrix300All[,1]<-c(0.42, 0.08, -0.07, -0.16, 0.984, 0.939, 0, 0.1,0.1,0.2,0.2,0.3,0.3,0.4,0.4)


Matrix300All[1,c(2,3,4)]<-Simulation(Iterations, SampleSize, "weib", 100, 2, 0, Matrix300All[1,1])
Matrix300All[2,c(2,3,4)]<-Simulation(Iterations, SampleSize, "weib", 100, 3, 0, Matrix300All[2,1])
Matrix300All[3,c(2,3,4)]<-Simulation(Iterations, SampleSize, "weib", 100, 4, 0, Matrix300All[3,1])
Matrix300All[4,c(2,3,4)]<-Simulation(Iterations, SampleSize, "weib", 100, 5, 0, Matrix300All[4,1])

Matrix300All[5,c(2,3,4)]<-Simulation(Iterations, SampleSize, "fnorm", 100, 0, 0, Matrix300All[5,1])
Matrix300All[6,c(2,3,4)]<-Simulation(Iterations, SampleSize, "lognorm", 100, 0, 0, Matrix300All[6,1])
Matrix300All[7,c(2,3,4)]<-Simulation(Iterations, SampleSize, "norm", 100, 0, 0, Matrix300All[7,1])
Matrix300All[8,c(2,3,4)]<-Simulation(Iterations, SampleSize, "normmixt", 100, 0.947, 0, Matrix300All[8,1])

Matrix300All[9,c(2,3,4)]<-Simulation(Iterations,  SampleSize, "normmixt", 100, 0.087, 0, Matrix300All[9,1])
Matrix300All[10,c(2,3,4)]<-Simulation(Iterations, SampleSize, "normmixt", 100, 0.885, 0, Matrix300All[10,1])
Matrix300All[11,c(2,3,4)]<-Simulation(Iterations, SampleSize, "normmixt", 100, 0.149, 0, Matrix300All[11,1])
Matrix300All[12,c(2,3,4)]<-Simulation(Iterations, SampleSize, "normmixt", 100, 0.809, 0, Matrix300All[12,1])

Matrix300All[13,c(2,3,4)]<-Simulation(Iterations, SampleSize, "normmixt", 100, 0.207, 0, Matrix300All[13,1])
Matrix300All[14,c(2,3,4)]<-Simulation(Iterations, SampleSize, "normmixt", 100, 0.708, 0, Matrix300All[14,1])
Matrix300All[15,c(2,3,4)]<-Simulation(Iterations, SampleSize, "normmixt", 100, 0.278, 0, Matrix300All[15,1])



Matrix300All

#module(envstats, unload=T)
#ooo<- MaxAsymCoef(400, "weib", 100, 2, 0)
#ooo